O Brasil destaca-se como o maior exportador de café no mundo e o segundo maior consumidor dessa bebida. Responsável por um terço da produção global, o Brasil mantém a liderança na produção de café há mais de 150 anos.
Entre as regiões produtoras, Minas Gerais é o maior estado produtor, contribuindo com cerca de 50% da produção nacional. Já São Paulo - conhecido pela sua tradição no cultivo de café, - produz exclusivamente Arábica.
Por ser uma paixão nacional, também somos o segundo maior consumidor no mundo, perdendo apenas para os EUA.
Em 2023 pela primeira vez na história, um produtor brasileiro, a São Mateus Agropecuária, conquistou o título “Best of the Best” no Prêmio Internacional Ernesto Illy. Esse reconhecimento é resultado de décadas de investimentos em técnicas sustentáveis e agricultura regenerativa. A fazenda foi pioneira em práticas como reuso da água e cuidado com o ecossistema, contribuindo para a qualidade do café brasileiro.
Fonte1: CONAB. Disponível em: https://www.abic.com.br/tudo-de-cafe/o-cafe-brasileiro-na-atualidade/. Acesso em: 02/12/23. Fonte2: Ansa. Disponível em: https://www.abic.com.br/tudo-de-cafe/o-cafe-brasileiro-na-atualidade/. Acesso em: 02/12/23. Fonte3: Disponível em: https://satocomunicacao.com.br/cafe-e-paixao-nacional-e-segunda-bebida-mais-consumida-pelos-brasileiros/. Acesso em: 02/12/23.
O objetivo principal deste projeto é desenvolver e validar um modelo de CLASSIFICÃO, capaz de distinguir cafés brasileiros de não-brasileiros com base em avaliações detalhadas fornecidas por Q Graders profissionais.
library(tidymodels)
library(tidyverse)
library(skimr)
library(ggplot2)
library(gridExtra)
library(patchwork)
library(reshape2)
library(car)
library(corrplot)
library(ranger)
library(MASS)
library(keras)
library(rsample)
library(yardstick)
library(FactoMineR)
library(factoextra)
library(doParallel)
library(caret)
library(pROC)
library(purrr)
library(gt)
library(ggrepel)
library(dplyr)
library(rmarkdown)
#2)Import da base
dados_originais <- read.csv("coffee_ratings_1.csv")
paged_table(dados_originais)
# Mostre o número de linhas e colunas do DataFrame
cat("Número de linhas:", nrow(dados_originais), "\n")
## Número de linhas: 1339
cat("Número de colunas:", ncol(dados_originais), "\n")
## Número de colunas: 43
#Verificando os dados missing
resultado <- dados_originais %>%
lapply(type_sum) %>%
as_tibble() %>%
pivot_longer(cols = 1:ncol(dados_originais),
names_to = "Coluna",
values_to = "Tipo") %>%
mutate(Contagem_NA = colSums(is.na(dados_originais))) %>%
arrange(desc(Contagem_NA))
resultado_df <- as.data.frame(resultado)
print(resultado_df)
## Coluna Tipo Contagem_NA
## 1 lot_number chr 1063
## 2 farm_name chr 359
## 3 mill chr 315
## 4 producer chr 231
## 5 altitude_low_meters dbl 230
## 6 altitude_high_meters dbl 230
## 7 altitude_mean_meters dbl 230
## 8 altitude chr 226
## 9 variety chr 226
## 10 color chr 218
## 11 company chr 209
## 12 processing_method chr 170
## 13 ico_number chr 151
## 14 region chr 59
## 15 harvest_year chr 47
## 16 owner chr 7
## 17 owner_1 chr 7
## 18 country_of_origin chr 1
## 19 quakers int 1
## 20 total_cup_points dbl 0
## 21 species chr 0
## 22 number_of_bags int 0
## 23 bag_weight chr 0
## 24 in_country_partner chr 0
## 25 grading_date chr 0
## 26 aroma dbl 0
## 27 flavor dbl 0
## 28 aftertaste dbl 0
## 29 acidity dbl 0
## 30 body dbl 0
## 31 balance dbl 0
## 32 uniformity dbl 0
## 33 clean_cup dbl 0
## 34 sweetness dbl 0
## 35 cupper_points dbl 0
## 36 moisture dbl 0
## 37 category_one_defects int 0
## 38 category_two_defects int 0
## 39 expiration chr 0
## 40 certification_body chr 0
## 41 certification_address chr 0
## 42 certification_contact chr 0
## 43 unit_of_measurement chr 0
# Dropando preditoras desinteressantes apos definiciao grupo
dados_tratados_1 <- dados_originais %>%
dplyr::select(- moisture,-cupper_points ,-category_two_defects ,-quakers ,-category_one_defects,-owner, -region, -farm_name, -lot_number, -mill, -ico_number, -company, -altitude, -producer, -number_of_bags, -bag_weight,
-in_country_partner, -harvest_year, -grading_date, -owner_1, -variety, -processing_method, -color,
-expiration, -certification_body, -certification_address, -certification_contact,-unit_of_measurement,-altitude_low_meters, -altitude_high_meters,-altitude_mean_meters)
#4.1) Dropando linhas missing da coluna region
dados_tratados_2 <- dados_tratados_1 %>%
dplyr::filter(country_of_origin != "Cote d?Ivoire")
# Excluindo o outlier em Honduras
dados_tratados_4 <- dados_tratados_2 %>%
filter(!(country_of_origin == "Honduras" & total_cup_points == 0))
#Estatistica
glimpse(dados_tratados_4)
## Rows: 1,336
## Columns: 12
## $ total_cup_points <dbl> 90.58, 89.92, 89.75, 89.00, 88.83, 88.83, 88.75, 88.…
## $ species <chr> "Arabica", "Arabica", "Arabica", "Arabica", "Arabica…
## $ country_of_origin <chr> "Ethiopia", "Ethiopia", "Guatemala", "Ethiopia", "Et…
## $ aroma <dbl> 8.67, 8.75, 8.42, 8.17, 8.25, 8.58, 8.42, 8.25, 8.67…
## $ flavor <dbl> 8.83, 8.67, 8.50, 8.58, 8.50, 8.42, 8.50, 8.33, 8.67…
## $ aftertaste <dbl> 8.67, 8.50, 8.42, 8.42, 8.25, 8.42, 8.33, 8.50, 8.58…
## $ acidity <dbl> 8.75, 8.58, 8.42, 8.42, 8.50, 8.50, 8.50, 8.42, 8.42…
## $ body <dbl> 8.50, 8.42, 8.33, 8.50, 8.42, 8.25, 8.25, 8.33, 8.33…
## $ balance <dbl> 8.42, 8.42, 8.42, 8.25, 8.33, 8.33, 8.25, 8.50, 8.42…
## $ uniformity <dbl> 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.…
## $ clean_cup <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, …
## $ sweetness <dbl> 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 9.3…
skim(dados_tratados_4)
| Name | dados_tratados_4 |
| Number of rows | 1336 |
| Number of columns | 12 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| numeric | 10 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| species | 0 | 1 | 7 | 7 | 0 | 2 | 0 |
| country_of_origin | 0 | 1 | 4 | 28 | 0 | 35 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| total_cup_points | 0 | 1 | 82.16 | 2.69 | 59.83 | 81.17 | 82.50 | 83.67 | 90.58 | ▁▁▁▇▁ |
| aroma | 0 | 1 | 7.57 | 0.32 | 5.08 | 7.42 | 7.58 | 7.75 | 8.75 | ▁▁▂▇▁ |
| flavor | 0 | 1 | 7.53 | 0.34 | 6.08 | 7.33 | 7.58 | 7.75 | 8.83 | ▁▂▇▃▁ |
| aftertaste | 0 | 1 | 7.41 | 0.35 | 6.17 | 7.25 | 7.42 | 7.60 | 8.67 | ▁▃▇▂▁ |
| acidity | 0 | 1 | 7.54 | 0.32 | 5.25 | 7.33 | 7.58 | 7.75 | 8.75 | ▁▁▃▇▁ |
| body | 0 | 1 | 7.52 | 0.31 | 5.08 | 7.33 | 7.50 | 7.67 | 8.58 | ▁▁▁▇▁ |
| balance | 0 | 1 | 7.52 | 0.35 | 5.25 | 7.33 | 7.50 | 7.75 | 8.75 | ▁▁▃▇▁ |
| uniformity | 0 | 1 | 9.84 | 0.49 | 6.00 | 10.00 | 10.00 | 10.00 | 10.00 | ▁▁▁▁▇ |
| clean_cup | 0 | 1 | 9.84 | 0.72 | 0.00 | 10.00 | 10.00 | 10.00 | 10.00 | ▁▁▁▁▇ |
| sweetness | 0 | 1 | 9.86 | 0.55 | 1.33 | 10.00 | 10.00 | 10.00 | 10.00 | ▁▁▁▁▇ |
#frequencia
table(dados_tratados_4$country_of_origin)
##
## Brazil Burundi
## 132 2
## China Colombia
## 16 183
## Costa Rica Ecuador
## 51 3
## El Salvador Ethiopia
## 21 44
## Guatemala Haiti
## 181 6
## Honduras India
## 52 14
## Indonesia Japan
## 20 1
## Kenya Laos
## 25 3
## Malawi Mauritius
## 11 1
## Mexico Myanmar
## 236 8
## Nicaragua Panama
## 26 4
## Papua New Guinea Peru
## 1 10
## Philippines Rwanda
## 5 1
## Taiwan Tanzania, United Republic Of
## 75 40
## Thailand Uganda
## 32 36
## United States United States (Hawaii)
## 10 73
## United States (Puerto Rico) Vietnam
## 4 8
## Zambia
## 1
table(dados_tratados_4$species)
##
## Arabica Robusta
## 1308 28
#Grande desbalanceamento de dados entre robusta e Arabica causando efeito de multicolinearidade
dados_tratados_4 <- dplyr::select(dados_tratados_4, -species)
# Balanceamento paises
# Tiblbe com a frequencia dos dados Paises
tabela_regioes <- dados_tratados_4 %>%
count(country_of_origin) %>%
rename(Contagem = n) %>%
group_by(country_of_origin) %>%
mutate(Percentual_Participacao = (Contagem / nrow(dados_tratados_4)) * 100) %>%
ungroup() %>%
arrange(desc(Contagem))
# Adicionar uma linha para somar os totais de contagem e participação
tabela_regioes <- tabela_regioes %>%
bind_rows(
tibble(
country_of_origin = "Total",
Contagem = sum(tabela_regioes$Contagem),
Percentual_Participacao = sum(tabela_regioes$Percentual_Participacao)
)
)
# Removendo a linha do 'Total' para o gráfico (se não quiser incluir o total)
tabela_regioes <- tabela_regioes %>% filter(country_of_origin != "Total")
# Reordenando os países com base na Contagem para exibição decrescente
tabela_regioes <- tabela_regioes %>%
mutate(country_of_origin = reorder(country_of_origin, -Contagem))
# Calculando o máximo da Contagem e Percentual_Participacao para a escala
max_contagem <- max(tabela_regioes$Contagem)
max_percentual <- max(tabela_regioes$Percentual_Participacao)
# Criando o gráfico
grafico <- ggplot(tabela_regioes, aes(x = country_of_origin, y = Contagem)) +
geom_bar(stat = "identity", fill = "blue") +
geom_line(aes(y = Percentual_Participacao * max_contagem / max_percentual, group = 1),
color = "red", linewidth = 0.5) +
geom_point(aes(y = Percentual_Participacao * max_contagem / max_percentual),
color = "red", size = 2) +
geom_text(aes(label = Contagem, y = Contagem), vjust = -0.5, size = 3.3) +
geom_text(aes(y = Percentual_Participacao * max_contagem / max_percentual, label = ifelse(Percentual_Participacao > 1, paste0(round(Percentual_Participacao, 0), "%"), "")),
color = "white", vjust = 2.5, size = 2.3) +
scale_y_continuous(sec.axis = sec_axis(~ . * max_percentual / max_contagem, name = "Percentual de Participação (%)")) +
labs(x = "Países", y = "Contagem") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Exibindo o gráfico
print(grafico)
#### 6.2) Nota Média por País
# 2.2) Nota media por pais
# Calculando a média de total_cup_points por país
media_pontuacao_por_pais <- dados_tratados_4 %>%
group_by(country_of_origin) %>%
summarise(media_pontuacao = mean(total_cup_points , na.rm = TRUE)) %>%
arrange(desc(media_pontuacao))
# Criando o gráfico com ggplot2
plot <- ggplot(media_pontuacao_por_pais, aes(x = reorder(country_of_origin, -media_pontuacao), y = media_pontuacao)) +
geom_bar(stat = "identity", fill = "steelblue") +
geom_text(aes(label = sprintf("%.2f", round(media_pontuacao, 2))),
position = position_dodge(width = 0.9), hjust = 1.1, size = 3.5, color = "white") +
theme_minimal() +
labs(title = "Média de Pontuação Total do Café por País", x = "País", y = "Média de Pontuação") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 0.5),
plot.title = element_text(size = 18, face = "bold"),
axis.title = element_text(size = 16)) +
coord_flip() # Usado para inverter as coordenadas e facilitar a leitura
# Plotando o gráfico
print(plot)
# Criar um único gráfico de dispersão para Total Cup Points vs paises
g <- ggplot(dados_tratados_4, aes(x = country_of_origin, y = total_cup_points)) +
geom_point(position = position_jitter(width = 0.2), alpha = 0.5) +
geom_smooth(aes(group = country_of_origin), method = "lm", se = FALSE, color = "blue") +
labs(title = "Dispersao do Total Cup Points por País") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Exibindo o gráfico
print(g)
## `geom_smooth()` using formula = 'y ~ x'
por Atributo
# 2.2) Nota media Brasil por atributo
# Filtrar os dados para incluir apenas as entradas do Brasil
dados_graf_brasil <- dados_tratados_4 %>%
filter(country_of_origin == "Brazil")
# Calcular a média de cada atributo para o Brasil
medias_brasil <- dados_graf_brasil %>%
summarise(
aroma = mean(aroma, na.rm = TRUE),
flavor = mean(flavor, na.rm = TRUE),
aftertaste = mean(aftertaste, na.rm = TRUE),
acidity = mean(acidity, na.rm = TRUE),
body = mean(body, na.rm = TRUE),
balance = mean(balance, na.rm = TRUE),
uniformity = mean(uniformity, na.rm = TRUE),
clean_cup = mean(clean_cup, na.rm = TRUE),
sweetness = mean(sweetness, na.rm = TRUE)
) %>%
mutate(country = "Brasil")
# Filtrar os dados para incluir apenas as entradas do Brasil
dados_graf_outros <- dados_tratados_4 %>%
filter(country_of_origin != "Brazil")
# Calcular a média de cada atributo para os outros países
medias_outros_paises <- dados_graf_outros %>%
summarise(
aroma = mean(aroma, na.rm = TRUE),
flavor = mean(flavor, na.rm = TRUE),
aftertaste = mean(aftertaste, na.rm = TRUE),
acidity = mean(acidity, na.rm = TRUE),
body = mean(body, na.rm = TRUE),
balance = mean(balance, na.rm = TRUE),
uniformity = mean(uniformity, na.rm = TRUE),
clean_cup = mean(clean_cup, na.rm = TRUE),
sweetness = mean(sweetness, na.rm = TRUE)
) %>%
mutate(country = "Outros Países")
# Unir os dados do Brasil com os dados dos outros países
dados_comparativos <- rbind(medias_brasil, medias_outros_paises) %>%
gather(attribute, mean_value, -country)
# Criando o gráfico com ggplot2
plot_comparativo <- ggplot(dados_comparativos, aes(x = attribute, y = mean_value, fill = country)) +
geom_col(position = position_dodge(), width = 0.8) +
geom_text(aes(label = sprintf("%.2f%%", mean_value)),
position = position_dodge(width = 0.8), vjust = 1.5, size = 3.7, angle = 45) +
theme_minimal() +
labs(title = "Comparação Brasil vs. Outros Países",
x = "Atributo", y = "Média") +
scale_fill_manual(values = c("Brasil" = "steelblue", "Outros Países" = "grey")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
plot.title = element_text(size = 18, face = "bold"),
axis.title = element_text(size = 16))
# Exibir o gráfico
print(plot_comparativo)
#### 6.5) Dipersão Total Notas vs Pais
# Criar um único gráfico de dispersão para Total Cup Points vs paises
g <- ggplot(dados_tratados_4, aes(x = country_of_origin, y = total_cup_points)) +
geom_point(position = position_jitter(width = 0.2), alpha = 0.5) +
geom_smooth(aes(group = country_of_origin), method = "lm", se = FALSE, color = "blue") +
labs(title = "Dispersao do Total Cup Points por País") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Exibindo o gráfico
print(g)
## `geom_smooth()` using formula = 'y ~ x'
dados_tratados_4 <- dplyr::select(dados_tratados_4, -total_cup_points)
#5) Checando as limpezas de dados antes e depois
dim_antes <- dim(dados_originais)
dim_depois <- dim(dados_tratados_4)
# Criando uma tibble para exibir as dimensões
dimensoes <- tibble(
"Checagem" = c("Antes", "Depois"),
"Linhas" = c(dim_antes[1], dim_depois[1]),
"Colunas" = c(dim_antes[2], dim_depois[2])
)
# Exibindo a tibble
print(dimensoes)
## # A tibble: 2 × 3
## Checagem Linhas Colunas
## <chr> <int> <int>
## 1 Antes 1339 43
## 2 Depois 1336 10
dados_final <- dados_tratados_4
# Exibindo a tibble
paged_table(dados_final)
# Definindo a semente para reprodutibilidade
set.seed(15)
Aqui precisamos usar uma técnica de Estratificacao para poder balancear os dados , ja que o número de observações para o Brasil estao desproporcionais comparado com os demais paises. Sem esse balanceamento o modelo ficará tendencioso, ou seja , ele será capaz de classificar com muita precisào o que nào e café brasileiro do que brasileiro , nao sendo esse nosso objetivo.
# Subamostragem dos dados de outros países para equilibrar com os dados do Brasil
dados_brazil <- dados_final %>% filter(country_of_origin == "Brazil")
dados_outros <- dados_final %>% filter(country_of_origin != "Brazil")
n_brazil <- nrow(dados_brazil)
dados_outros_subamostrados <- dados_outros %>%
slice_sample(n = n_brazil)
# Combinando os dados subamostrados com os dados do Brasil
dados_combinados <- bind_rows(dados_brazil, dados_outros_subamostrados)
# Redefinindo a divisão de treino e teste com os dados balanceados
split_combinado <- initial_split(dados_combinados, prop = 0.7)
treinamento <- training(split_combinado)
teste <- testing(split_combinado)
# Checando o Balanceamento____________________________________________________
n_brasil_treinamento <- treinamento %>%
filter(country_of_origin == "Brazil") %>%
summarize(count = n()) %>%
pull(count)
# Contando os registros do Brasil no conjunto de teste
n_brasil_teste <- teste %>%
filter(country_of_origin == "Brazil") %>%
summarize(count = n()) %>%
pull(count)
# Criando uma tibble para mostrar os resultados
resumo_brasil <- tibble(
Conjunto = c("Treinamento", "Teste"),
Total_Brasil = c(n_brasil_treinamento, n_brasil_teste)
)
# Exibindo o resumo final
print(resumo_brasil)
## # A tibble: 2 × 2
## Conjunto Total_Brasil
## <chr> <int>
## 1 Treinamento 95
## 2 Teste 37
print(nrow(treinamento))
## [1] 184
print(nrow(teste))
## [1] 80
# Checando o Balanceamento____________________________________________________
# Criar a nova variável alvo para classificação binária como fator
treinamento$from_brazil <- as.factor(ifelse(treinamento$country_of_origin == "Brazil", 1, 0))
teste$from_brazil <- as.factor(ifelse(teste$country_of_origin == "Brazil", 1, 0))
# Remover a coluna original 'country_of_origin'
treinamento$country_of_origin <- NULL
teste$country_of_origin <- NULL
# Preparando a receita com a nova variável alvo
receita <- recipe(from_brazil ~ ., data = treinamento) %>%
step_normalize(all_numeric())
receita_prep <- prep(receita)
treinamento_proc <- bake(receita_prep, new_data = NULL)
teste_proc <- bake(receita_prep, new_data = teste)
fit_glm <- logistic_reg() %>%
set_engine("glm") %>%
set_mode("classification")
logistica_simples <- fit(fit_glm, from_brazil ~ ., data = treinamento_proc)
# Fazendo previsões no conjunto de teste e combinando com os dados de teste
logistica <- predict(logistica_simples, teste_proc, type = "class") %>%
bind_cols(teste_proc) %>%
rename(Predito = .pred_class)
# Verificando se a coluna 'from_brazil' está presente
if (!"from_brazil" %in% names(logistica)) {
logistica$from_brazil <- teste$from_brazil
}
# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
logistica <- logistica %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))
#Definindo o modelo de regressão logística com penalidade Ridge
ridge <- logistic_reg(penalty = tune(), mixture = 1) %>%
set_engine("glmnet") %>%
set_mode("classification")
# Validação cruzada para ajuste do hiperparâmetro em 10 lotes
ridge_cv_split <- vfold_cv(treinamento_proc, v = 10)
doParallel::registerDoParallel() # Paralelizando os próximos comandos
# Definição do espaço de busca para o hiperparâmetro 'penalty'
val_gride_ridge <- grid_regular(penalty(range = c(0.001, 1)), levels = 20)
# Ajuste do hiperparâmetro usando validação cruzada
ridge_lambda_tune <- tune_grid(
ridge,
receita,
resamples = ridge_cv_split,
grid = val_gride_ridge,
metrics = metric_set(accuracy),
control = control_grid(save_pred = TRUE)
)
# Selecionando a melhor combinação de hiperparâmetros
ridge_best <- ridge_lambda_tune %>%
select_best("accuracy")
# Finalizando o modelo com os melhores hiperparâmetros
fit_ridge <- finalize_model(ridge, parameters = ridge_best) %>%
fit(from_brazil ~ ., data = treinamento_proc)
# Previsões no conjunto de teste
ridge_fit_tunado <- predict(fit_ridge, teste_proc, type = "class") %>%
bind_cols(teste_proc) %>%
rename(Predito = .pred_class)
# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
ridge_fit_tunado <- ridge_fit_tunado %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))
# Definindo o modelo de regressão logística com penalidade Ridge
lasso <- logistic_reg(penalty = tune(), mixture = 0) %>%
set_engine("glmnet") %>%
set_mode("classification")
# Validação cruzada para ajuste do hiperparâmetro em 10 lotes
ridge_cv_split <- vfold_cv(treinamento_proc, v = 10)
doParallel::registerDoParallel() # Paralelizando os próximos comandos
# Definição do espaço de busca para o hiperparâmetro 'penalty'
val_gride_lasso <- grid_regular(penalty(range = c(0.001, 1)), levels = 20)
# Ajuste do hiperparâmetro usando validação cruzada
lasso_lambda_tune <- tune_grid(
lasso,
receita,
resamples = ridge_cv_split,
grid = val_gride_lasso,
metrics = metric_set(accuracy),
control = control_grid(save_pred = TRUE)
)
# Selecionando a melhor combinação de hiperparâmetros
lasso_best <- lasso_lambda_tune %>%
select_best("accuracy")
# Finalizando o modelo com os melhores hiperparâmetros
fit_lasso <- finalize_model(lasso, parameters = lasso_best) %>%
fit(from_brazil ~ ., data = treinamento_proc)
# Previsões no conjunto de teste
lasso_fit_tunado <- predict(fit_lasso, teste_proc, type = "class") %>%
bind_cols(teste_proc) %>%
rename(Predito = .pred_class)
# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
lasso_fit_tunado <- lasso_fit_tunado %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))
# Definindo o modelo de árvore de decisão
arvore <- decision_tree(tree_depth = tune(), cost_complexity = tune()) %>%
set_engine("rpart") %>%
set_mode("classification")
# Paralelizando os próximos comandos
doParallel::registerDoParallel()
# Definição do espaço de busca para os hiperparâmetros 'tree_depth' e 'cost_complexity'
val_grid_arvore <- grid_regular(
tree_depth(range = c(1, 10)),
cost_complexity(range = c(0.001, 0.1)),
levels = 10
)
# Criando o objeto de validação cruzada
cv_split <- vfold_cv(treinamento_proc, v = 10)
# Ajuste do hiperparâmetro usando validação cruzada
arvore_tune <- tune_grid(
arvore,
receita,
resamples = cv_split,
grid = val_grid_arvore,
metrics = metric_set(accuracy)
)
# Coletando as métricas do modelo ajustado
arvore_metrics <- arvore_tune %>%
collect_metrics()
# Selecionando a melhor combinação de hiperparâmetros com base na acurácia
arvore_best <- arvore_tune %>%
select_best("accuracy")
# Finalizando o modelo com os melhores hiperparâmetros
fit_arvore_final <- finalize_model(arvore, parameters = arvore_best) %>%
fit(from_brazil ~ ., data = treinamento_proc)
# Previsões no conjunto de teste
arvore_fit_tunado <- predict(fit_arvore_final, teste_proc, type = "class") %>%
bind_cols(teste_proc) %>%
rename(Predito = .pred_class)
# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
arvore_fit_tunado <- arvore_fit_tunado %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))
# Define o modelo XGBoost com parâmetros para serem tunados
xgboost_model <- boost_tree(
trees = tune(), # Número de árvores
tree_depth = tune(), # Profundidade da árvore
min_n = tune(), # Número mínimo de observações
loss_reduction = tune(), # Redução de perda necessária para fazer uma divisão adicional
mtry = tune(), # Número de variáveis a considerar em cada divisão
learn_rate = tune() # Taxa de aprendizado
) %>%
set_engine("xgboost") %>%
set_mode("classification")
# Criando o objeto de validação cruzada
cv_split_boost <- vfold_cv(treinamento_proc, v = 10)
# Executa a sintonia (tuning) do modelo
xgboost_tune <- tune_grid(
xgboost_model,
receita,
resamples = cv_split_boost,
grid = 30,
metrics = metric_set(accuracy)
)
## i Creating pre-processing data to finalize unknown parameter: mtry
# Coleta as métricas calculadas
xgboost_metrics <- xgboost_tune %>%
collect_metrics()
# Seleciona a melhor combinação de hiperparâmetros
xgboost_best <- xgboost_tune %>%
select_best("accuracy")
# Finaliza o modelo com os melhores hiperparâmetros
fit_xgboost <- finalize_model(xgboost_model, parameters = xgboost_best) %>%
fit(from_brazil ~ ., data = treinamento_proc)
# Faz previsões no conjunto de teste
xgboost_fit_tunado <- predict(fit_xgboost, teste_proc, type = "class") %>%
bind_cols(teste_proc) %>%
rename(Predito = .pred_class) %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == "1", "Brasileiro", "Não Brasileiro"))
# Preparação dos dados para a Rede Neural **************************************
# Seleciona as variáveis preditoras e converte para matriz
X_trn <- treinamento_proc %>%
dplyr::select(-from_brazil) %>%
as.matrix()
X_tst <- teste_proc %>%
dplyr::select(-from_brazil) %>%
as.matrix()
# Normalização dos dados
X_trn <- scale(X_trn)
X_tst <- scale(X_tst, center = attr(X_trn, "scaled:center"), scale = attr(X_trn, "scaled:scale"))
# Converter as classes de destino para one-hot encoding
num_classes <- length(unique(treinamento_proc$from_brazil))
Y_trn <- to_categorical(as.numeric(factor(treinamento_proc$from_brazil)) - 1, num_classes)
Y_tst <- to_categorical(as.numeric(factor(teste_proc$from_brazil)) - 1, num_classes)
# Construção da Rede Neural ****************************************************
net <- keras_model_sequential() %>%
layer_dense(units = 64, activation = "relu", input_shape = ncol(X_trn)) %>%
layer_dense(units = 32, activation = "relu") %>%
layer_dense(units = 16, activation = "relu") %>%
layer_dense(units = num_classes, activation = "sigmoid")
net <- compile(
net,
loss = "binary_crossentropy",
optimizer = "adam",
metrics = "accuracy"
)
summary(net)
## Model: "sequential"
## ________________________________________________________________________________
## Layer (type) Output Shape Param #
## ================================================================================
## dense_3 (Dense) (None, 64) 640
## dense_2 (Dense) (None, 32) 2080
## dense_1 (Dense) (None, 16) 528
## dense (Dense) (None, 2) 34
## ================================================================================
## Total params: 3282 (12.82 KB)
## Trainable params: 3282 (12.82 KB)
## Non-trainable params: 0 (0.00 Byte)
## ________________________________________________________________________________
# Treinamento da Rede Neural ****************************************************
history <- fit(net, X_trn, Y_trn, batch_size = 50, epochs = 30, validation_split = 0.2)
## Epoch 1/30
## 3/3 - 1s - loss: 0.7663 - accuracy: 0.4150 - val_loss: 0.8425 - val_accuracy: 0.4054 - 603ms/epoch - 201ms/step
## Epoch 2/30
## 3/3 - 0s - loss: 0.7310 - accuracy: 0.4898 - val_loss: 0.8241 - val_accuracy: 0.4324 - 49ms/epoch - 16ms/step
## Epoch 3/30
## 3/3 - 0s - loss: 0.7093 - accuracy: 0.5510 - val_loss: 0.8081 - val_accuracy: 0.5405 - 37ms/epoch - 12ms/step
## Epoch 4/30
## 3/3 - 0s - loss: 0.6980 - accuracy: 0.5646 - val_loss: 0.7919 - val_accuracy: 0.5405 - 36ms/epoch - 12ms/step
## Epoch 5/30
## 3/3 - 0s - loss: 0.6879 - accuracy: 0.5714 - val_loss: 0.7754 - val_accuracy: 0.5676 - 36ms/epoch - 12ms/step
## Epoch 6/30
## 3/3 - 0s - loss: 0.6790 - accuracy: 0.5918 - val_loss: 0.7614 - val_accuracy: 0.5135 - 36ms/epoch - 12ms/step
## Epoch 7/30
## 3/3 - 0s - loss: 0.6711 - accuracy: 0.6190 - val_loss: 0.7493 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 8/30
## 3/3 - 0s - loss: 0.6640 - accuracy: 0.6122 - val_loss: 0.7419 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 9/30
## 3/3 - 0s - loss: 0.6587 - accuracy: 0.6122 - val_loss: 0.7355 - val_accuracy: 0.4865 - 36ms/epoch - 12ms/step
## Epoch 10/30
## 3/3 - 0s - loss: 0.6532 - accuracy: 0.6259 - val_loss: 0.7338 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 11/30
## 3/3 - 0s - loss: 0.6486 - accuracy: 0.6395 - val_loss: 0.7335 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 12/30
## 3/3 - 0s - loss: 0.6446 - accuracy: 0.6395 - val_loss: 0.7306 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 13/30
## 3/3 - 0s - loss: 0.6416 - accuracy: 0.6395 - val_loss: 0.7322 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 14/30
## 3/3 - 0s - loss: 0.6392 - accuracy: 0.6531 - val_loss: 0.7338 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 15/30
## 3/3 - 0s - loss: 0.6362 - accuracy: 0.6599 - val_loss: 0.7352 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 16/30
## 3/3 - 0s - loss: 0.6344 - accuracy: 0.6667 - val_loss: 0.7365 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 17/30
## 3/3 - 0s - loss: 0.6321 - accuracy: 0.6599 - val_loss: 0.7387 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 18/30
## 3/3 - 0s - loss: 0.6306 - accuracy: 0.6667 - val_loss: 0.7434 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 19/30
## 3/3 - 0s - loss: 0.6282 - accuracy: 0.6599 - val_loss: 0.7497 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 20/30
## 3/3 - 0s - loss: 0.6259 - accuracy: 0.6667 - val_loss: 0.7515 - val_accuracy: 0.5135 - 35ms/epoch - 12ms/step
## Epoch 21/30
## 3/3 - 0s - loss: 0.6242 - accuracy: 0.6667 - val_loss: 0.7589 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
## Epoch 22/30
## 3/3 - 0s - loss: 0.6232 - accuracy: 0.6735 - val_loss: 0.7651 - val_accuracy: 0.4595 - 34ms/epoch - 11ms/step
## Epoch 23/30
## 3/3 - 0s - loss: 0.6214 - accuracy: 0.6667 - val_loss: 0.7712 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 24/30
## 3/3 - 0s - loss: 0.6201 - accuracy: 0.6803 - val_loss: 0.7754 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 25/30
## 3/3 - 0s - loss: 0.6202 - accuracy: 0.6803 - val_loss: 0.7769 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 26/30
## 3/3 - 0s - loss: 0.6187 - accuracy: 0.6735 - val_loss: 0.7850 - val_accuracy: 0.4595 - 34ms/epoch - 11ms/step
## Epoch 27/30
## 3/3 - 0s - loss: 0.6173 - accuracy: 0.6735 - val_loss: 0.7917 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 28/30
## 3/3 - 0s - loss: 0.6166 - accuracy: 0.6735 - val_loss: 0.8002 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 29/30
## 3/3 - 0s - loss: 0.6177 - accuracy: 0.6735 - val_loss: 0.8054 - val_accuracy: 0.4595 - 35ms/epoch - 12ms/step
## Epoch 30/30
## 3/3 - 0s - loss: 0.6165 - accuracy: 0.6735 - val_loss: 0.8137 - val_accuracy: 0.4865 - 35ms/epoch - 12ms/step
# Previsões e Avaliação *********************************************************
y_hat_net <- predict(net, X_tst)
## 3/3 - 0s - 58ms/epoch - 19ms/step
predicted_classes <- apply(y_hat_net, 1, which.max) - 1 # Obter as classes previstas
# Criando um dataframe temporário para as previsões da Rede Neural
previsoes_rede_neural <- cbind(teste_proc, Predito = as.factor(predicted_classes))
# Convertendo '0' e '1' para 'Não Brasileiro' e 'Brasileiro'
previsoes_rede_neural <- previsoes_rede_neural %>%
mutate(Observado = ifelse(from_brazil == "1", "Brasileiro", "Não Brasileiro"),
Predito = ifelse(Predito == 1, "Brasileiro", "Não Brasileiro"))
# 10) Matriz de Confusão
criar_df_confusao <- function(predicoes, observado, nome_modelo) {
# Definindo os níveis na ordem desejada com "Brasileiro" primeiro
niveis <- c("Brasileiro", "Não Brasileiro")
# Convertendo para fatores com os níveis definidos
predicoes <- factor(predicoes, levels = niveis)
observado <- factor(observado, levels = niveis)
cm <- confusionMatrix(predicoes, observado)
df <- cbind(melt(as.matrix(cm$table)), Modelo = nome_modelo)
return(df)
}
# Convertendo previsões e observações para fatores com os mesmos níveis
logistica$Predito <- factor(logistica$Predito)
logistica$Observado <- factor(logistica$Observado)
df_logistica <- criar_df_confusao(logistica$Predito, logistica$Observado, "Logística")
ridge_fit_tunado$Predito <- factor(ridge_fit_tunado$Predito)
ridge_fit_tunado$Observado <- factor(ridge_fit_tunado$Observado)
df_ridge <- criar_df_confusao(ridge_fit_tunado$Predito, ridge_fit_tunado$Observado, "Ridge Tunado")
lasso_fit_tunado$Predito <- factor(lasso_fit_tunado$Predito)
lasso_fit_tunado$Observado <- factor(lasso_fit_tunado$Observado)
df_lasso <- criar_df_confusao(lasso_fit_tunado$Predito, lasso_fit_tunado$Observado, "Lasso Tunado")
arvore_fit_tunado$Predito <- factor(arvore_fit_tunado$Predito)
arvore_fit_tunado$Observado <- factor(arvore_fit_tunado$Observado)
df_arvore <- criar_df_confusao(arvore_fit_tunado$Predito, arvore_fit_tunado$Observado, "Árvore de Decisão")
xgboost_fit_tunado$Predito <- factor(xgboost_fit_tunado$Predito)
xgboost_fit_tunado$Observado <- factor(xgboost_fit_tunado$Observado)
df_xgboost <- criar_df_confusao(xgboost_fit_tunado$Predito, xgboost_fit_tunado$Observado, "XGBoost")
# Para a Rede Neural, ajustando as previsões
predicted_classes_net <- ifelse(previsoes_rede_neural$Predito == 1, "Brasileiro", "Não Brasileiro")
predicted_classes_net <- factor(predicted_classes_net)
previsoes_rede_neural$Observado <- factor(previsoes_rede_neural$Observado)
df_rede_neural <- criar_df_confusao(predicted_classes_net, previsoes_rede_neural$Observado, "Rede Neural")
# Combinando todos os dataframes
df_all <- rbind(df_logistica, df_ridge, df_lasso, df_arvore, df_xgboost, df_rede_neural)
ggplot(data = df_all, aes(x = Prediction, y = Reference)) +
geom_tile(aes(fill = value), color = 'darkblue') +
geom_text(aes(label = sprintf("%d", value)), vjust = 1, size = 3) +
scale_fill_gradient(low = "lightblue", high = "steelblue") +
scale_y_discrete(limits = rev(levels(df_all$Reference))) + # Isso irá inverter a ordem no eixo Y
theme_minimal() +
labs(title = "Comparação de Matrizes de Confusão", x = "Predição", y = "Real", fill = "Contagem") +
theme(axis.text.x = element_text(angle = 45, hjust = 1), axis.text.y = element_text()) +
facet_wrap(~ Modelo, ncol = 3)
# Para o modelo Logistica
logistica_probs <- predict(logistica_simples, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo Ridge
ridge_probs <- predict(fit_ridge, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo Lasso
lasso_probs <- predict(fit_lasso, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo de Árvore de Decisão
arvore_probs <- predict(fit_arvore_final, teste_proc, type = "prob")[[".pred_1"]]
# Para o modelo XGBoost
xgboost_probs <- predict(fit_xgboost, teste_proc, type = "prob")[[".pred_1"]]
# Para a Rede Neural
rede_neural_probs <- predict(net, X_tst)[, 1]
## 3/3 - 0s - 12ms/epoch - 4ms/step
resultados_preds <- bind_rows(
tibble(
modelo = "Logística",
marcador = logistica_probs,
resposta = as.factor(teste$from_brazil)
),
tibble(
modelo = "Ridge",
marcador = ridge_probs,
resposta = as.factor(teste$from_brazil)
),
tibble(
modelo = "Lasso",
marcador = lasso_probs,
resposta = as.factor(teste$from_brazil)
),
tibble(
modelo = "Árvore de Decisão",
marcador = arvore_probs,
resposta = as.factor(teste$from_brazil)
),
tibble(
modelo = "XGBoost",
marcador = xgboost_probs,
resposta = as.factor(teste$from_brazil)
),
tibble(
modelo = "Neural",
marcador = rede_neural_probs,
resposta = as.factor(teste$from_brazil)
)
)
# Calculando a curva ROC para cada modelo e plotando
roc_plot <- resultados_preds %>%
group_by(modelo) %>%
roc_curve(resposta, marcador, event_level = "second") %>%
autoplot() +
theme(legend.position = "right")
print(roc_plot)
( Acurácia, Precisão , Recall e F1 - Score)
#Avaliando metricas _
lista_modelos <- list(df_logistica, df_ridge, df_lasso, df_arvore, df_xgboost, df_rede_neural)
# Calculando as métricas diretamente dos dados fornecidos
metricas_modelos <- lapply(lista_modelos, function(df) {
verdadeiro_positivo <- df$value[df$Prediction == "Brasileiro" & df$Reference == "Brasileiro"]
falso_positivo <- df$value[df$Prediction == "Não Brasileiro" & df$Reference == "Brasileiro"]
falso_negativo <- df$value[df$Prediction == "Brasileiro" & df$Reference == "Não Brasileiro"]
verdadeiro_negativo <- df$value[df$Prediction == "Não Brasileiro" & df$Reference == "Não Brasileiro"]
acuracia <- (verdadeiro_positivo + verdadeiro_negativo) / sum(df$value)
precisao <- verdadeiro_positivo / (verdadeiro_positivo + falso_positivo)
recall <- verdadeiro_positivo / (verdadeiro_positivo + falso_negativo)
f1 <- ifelse((precisao + recall) == 0, 0, 2 * (precisao * recall) / (precisao + recall))
tibble(
Modelo = unique(df$Modelo),
Acuracia = acuracia,
Precisao = precisao,
Recall = recall,
F1 = f1
)
}) %>% bind_rows()
# Tratando NA, convertendo métricas para porcentagem e ordenando por acurácia
metricas_modelos <- metricas_modelos %>%
mutate(
across(c(Acuracia, Precisao, Recall, F1), ~ifelse(is.na(.), 0, round(. * 100, 0)))
) %>%
arrange(desc(Acuracia))
# Formatação dos dados para a visualização com gt
gt_table <- metricas_modelos %>%
gt() %>%
tab_header(
title = "Métricas de Desempenho dos Modelos"
) %>%
cols_label(
Modelo = "Modelo",
Acuracia = "Acurácia",
Precisao = "Precisão",
Recall = "Recall",
F1 = "F1 Score"
) %>%
fmt_number(
columns = c("Acuracia", "Precisao", "Recall", "F1"),
decimals = 0,
pattern = "{x}%"
) %>%
tab_options(
heading.background.color = "gray",
column_labels.font.size = 12,
data_row.padding = px(10)
) %>%
tab_style(
style = cell_fill(color = "lightblue"),
locations = cells_body(
columns = c("Acuracia", "Precisao", "Recall", "F1")
)
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels(
columns = c("Modelo", "Acuracia", "Precisao", "Recall", "F1")
)
)
# Exibindo a tabela
gt_table
| Métricas de Desempenho dos Modelos | ||||
| Modelo | Acurácia | Precisão | Recall | F1 Score |
|---|---|---|---|---|
| XGBoost | 65% | 65% | 62% | 63% |
| Rede Neural | 54% | 0% | 0% | 0% |
| Logística | 50% | 62% | 47% | 53% |
| Lasso Tunado | 48% | 81% | 46% | 59% |
| Ridge Tunado | 46% | 100% | 46% | 63% |
| Árvore de Decisão | 46% | 100% | 46% | 63% |
A ánalise foi feita para os países mais avaliados em número de observacões dentro das observacões ( Top 10)
# Passando os dados tratados para o nova variavel
dados_clust <- dados_final
# Contando os registros para cada país para levar somente o top 10 mais avaliados em frequencia
country_counts <- dados_clust %>%
count(country_of_origin) %>%
arrange(desc(n))
# Selecionando os top 6 países com mais registros
top_countries <- head(country_counts, 10)
# Filtrando 'dados_clust' para incluir apenas os top 6 países
filtered_clust <- dados_clust %>%
filter(country_of_origin %in% top_countries$country_of_origin)
# Agrupando e somando os atributos para os países filtrados
df <- filtered_clust %>%
group_by(country_of_origin) %>%
summarise(across(c(aroma, flavor, aftertaste, acidity, body, balance, uniformity, clean_cup, sweetness), ~sum(., na.rm = TRUE)))
# Definir a coluna 'regiao' como row names e remover a coluna do dataframe
df <- df %>%
column_to_rownames(var = "country_of_origin")
X <- scale(df,
center = TRUE, # centraliza os dados
scale = TRUE) # escalona os dados (pois estao em medidas diferentes)
pca <- prcomp(X) # aplica o PCA
# Todos os atributos para PC1 tem o mesmo peso nao mostrando atributos que se sobressaem
pca$rotation <- -pca$rotation # troca o sinal das cargas
pca$x <- -pca$x # troca o sinal dos scores
Phi <- pca$rotation # matriz de cargas
head(Phi)
## PC1 PC2 PC3 PC4 PC5 PC6
## aroma 0.3333679 0.01157345 0.08366266 -0.04978241 0.4417232 -0.15801712
## flavor 0.3333591 0.16693211 0.16272529 -0.05493036 -0.6843610 0.14748734
## aftertaste 0.3333299 0.34140572 -0.37119469 -0.36685208 -0.1632435 -0.26629708
## acidity 0.3333312 0.01458445 0.79592318 0.14359109 0.1305645 0.01362871
## body 0.3333602 0.18745743 0.10937622 -0.12224716 -0.1729299 -0.06988413
## balance 0.3333115 0.50462916 -0.23607065 0.12983608 0.4671684 0.30713836
## PC7 PC8 PC9
## aroma -0.17416192 0.61170571 0.50462770
## flavor -0.08661671 0.51395144 -0.26411654
## aftertaste -0.47601668 -0.38692467 0.16156189
## acidity -0.26335184 -0.37533772 -0.08426192
## body 0.77360425 -0.20443827 0.38958452
## balance 0.16453901 0.02848641 -0.47032037
Z <- pca$x # matriz de scores
head(Z)
## PC1 PC2 PC3 PC4 PC5
## Brazil 1.096309 0.007548140 -0.011464123 -0.004428219 -0.0152341163
## Colombia 3.347607 0.057000370 -0.031358419 0.007881091 0.0054410453
## Costa Rica -2.357759 -0.003584463 -0.003792201 -0.003037407 0.0083292517
## Ethiopia -2.594418 0.031025358 0.019949117 -0.006628561 -0.0043695017
## Guatemala 3.149011 -0.006995421 0.027747603 0.038912421 0.0005027964
## Honduras -2.365499 -0.033384541 -0.003929353 0.006117134 -0.0008986453
## PC6 PC7 PC8 PC9
## Brazil 0.002174996 0.0001237894 -0.0005462686 -2.258589e-05
## Colombia 0.001064819 0.0002477694 0.0001204255 1.233801e-05
## Costa Rica 0.001235835 -0.0035754693 0.0012226895 -3.276435e-05
## Ethiopia 0.001609311 -0.0046573502 0.0003097058 2.670071e-05
## Guatemala -0.000535794 -0.0001326772 -0.0001578867 -5.862193e-06
## Honduras 0.003455312 0.0049293446 0.0018220648 1.308488e-05
#Analises PCA_______________-
# Resumo do PCA para verificar a variação explicada por cada componente principal
summary(pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.9997 0.03520 0.01884 0.01808 0.006653 0.003218
## Proportion of Variance 0.9998 0.00014 0.00004 0.00004 0.000000 0.000000
## Cumulative Proportion 0.9998 0.99992 0.99996 0.99999 1.000000 1.000000
## PC7 PC8 PC9
## Standard deviation 0.003044 0.00126 1.786e-05
## Proportion of Variance 0.000000 0.00000 0.000e+00
## Cumulative Proportion 1.000000 1.00000 1.000e+00
# Obtendo as cargas para cada atributo nos componentes principais
loadings <- abs(pca$rotation)
# Para identificar os atributos mais importantes para o primeiro componente principal (PC1)
atributos_importantes_pc1 <- sort(loadings[,1], decreasing = TRUE)
# Imprimindo os atributos mais importantes para PC1
print(atributos_importantes_pc1)
## aroma body flavor clean_cup acidity uniformity aftertaste
## 0.3333679 0.3333602 0.3333591 0.3333313 0.3333312 0.3333308 0.3333299
## balance sweetness
## 0.3333115 0.3332782
# Convertendo a matriz Phi em uma tibble, preservando os nomes das linhas
contribuicoes <- as_tibble(Phi, rownames = "country_of_origin") %>%
mutate(Contribuicao_PC1 = (PC1^2) * 100 / sum(PC1^2)) %>%
dplyr::select(country_of_origin, Contribuicao_PC1)
# Visualizando as contribuições
print(contribuicoes)
## # A tibble: 9 × 2
## country_of_origin Contribuicao_PC1
## <chr> <dbl>
## 1 aroma 11.1
## 2 flavor 11.1
## 3 aftertaste 11.1
## 4 acidity 11.1
## 5 body 11.1
## 6 balance 11.1
## 7 uniformity 11.1
## 8 clean_cup 11.1
## 9 sweetness 11.1
# o grafico abaixo mostra o percentual explicado da variancia de cada componente
fviz_eig(pca, addlabels = TRUE) +
labs(x = "Componente Principal",
y = "Percentual explicado da variância")
# Grafico mostra a PCA explicando quase toda a variancia da base
# biplot do factoextra
fviz_pca_biplot(pca,
repel = TRUE,
xlab = "PC1 - Intensidade do Aroma",
ylab = "PC2 - Qualidade do Sabor",
labelsize = 3,
geom = c("point", "text"),
addEllipses = FALSE) +
theme(plot.margin = margin(0, 0, 0, 0, "cm"),
text = element_text(size = 10)) +
theme(axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 12),
axis.text.x = element_text(size = 8),
axis.text.y = element_text(size = 8))
Brasil: Está próximo do centro no eixo PC1, indicando uma intensidade de aroma moderada. Os atributos “body”, “balance” e “aftertaste” apontam na direção do Brasil, sugerindo que esses atributos são fortes no café brasileiro.
Nossa expectativa era alcançar uma precisão superior a 65% na identificação do café brasileiro. Contudo, enfrentamos desafios significativos devido ao desbalanceamento dos dados. Esse obstáculo nos proporcionou um valioso aprendizado sobre como manejar eficientemente tais problemas. Empregamos e descartamos múltiplos modelos de regressão altamente colinearizados, e adotamos diversas abordagens e técnicas para mitigar essas dificuldades. Essa jornada nos levou à conclusão de que uma acurácia de 65% e a utilização do método XGBoost representam a melhor performance possível para determinar se um café é brasileiro ou não.